clear all;close all;clc
fs = 44.1e3;
f1 = 200;
f2 = 1000;
f3 = 12000;

One1_Two2 = 1; % plot binaural HRTFs
% One1_Two2 = 2; % 2 plot dual HRTFs

folder = 'standard_hrir_database/';
SN = {'003';'008';'009';'010';'011';'012';'015';'017';'018';'019';'020';'021';'027';...
    '028';'033';'040';'044';'048';'051';'058';'059';'060';'061';'065';'119';'124';...
    '126';'127';'131';'133';'134';'135';'137';'147';'148';'152';'153';'154';'155';'156'};
X_ILD = [];X_ITD = [];
d = [];
c_row = 0;


xv = 2:2:24; % horizontal 1-25 locations ; official (12)
% xv = [1 8 15];
yv = 1:5:50; % vertical 1-50 locations; official
% yv = [4 10 26];
x_loc = [-80:15:-65 -55 -45:5:45 55 65:15:80]; %Azimuths
y_loc = -45:5.625:230.625; %Elevations

length(xv)
length(yv)
NC = length(xv)*length(yv)
c = 0;
for azn=xv
    c_row=c_row+1;
    c_col = 0;
    
    for eln=yv
        c_col=c_col+1;
        c = c+1;
        if One1_Two2 == 1
            Binaural_vectors{c} = [];
            Binaural_vectors{c} = [];
        else
            Binaural_vectors{c_row,c_col} = [];
            Binaural_vectors{c_row,c_col} = [];
            Monaural_vectors{c_row,c_col} = [];
            Monaural_vectors{c_row,c_col} = [];
        end
    end
end

for isub = 1:length(SN)
    isub
    c = 0;
    c_row = 0;
    eval(['load ' folder 'subject_' SN{isub} '/hrir_final'])
    for azn=xv
        c_col = 0;
        c_row=c_row+1;
        temp_vect = [];
        temp_mon = [];
        for eln=yv
            c_col=c_col+1;
            h_l = squeeze(hrir_l(azn,eln,:));
            h_r = squeeze(hrir_r(azn,eln,:));
            Nfft = length(h_l);
            fsig = (1:Nfft)./Nfft.*fs;
            N1 = find(fsig>f1);N1 = N1(1);
            N2 = find(fsig<=f2);N2 = N2(end);
            N3 = find(fsig<=f3);N3 = N3(end);
            yl = abs(fft(h_l,Nfft))';
            yr = abs(fft(h_r,Nfft))';
            phasel = angle(fft(h_l,Nfft))';
            phaser = angle(fft(h_r,Nfft))';
            
            if One1_Two2==1
                phrtf =  [(phaser(N1:N2)-phasel(N1:N2)) yr(N2:N3)-yl(N2:N3)];
                temp_vect=[temp_vect
                    phrtf];
            else        
                phrtf =  [(phaser(N1:N2)-phasel(N1:N2)) yr(N2:N3)-yl(N2:N3)];
                temp_vect=[temp_vect
                    phrtf];
                temp_mon=[temp_mon
                    yl(N1:N3)];
            end
        end % elevation finished
        if One1_Two2==1 % do the correction
           temp_vect=f_HRTF_wrap(temp_vect);        
            m = 0;
            for kk = c+1:c+length(yv)
                m = m+1;
                Binaural_vectors{kk}=[Binaural_vectors{kk}
                    temp_vect(m,:)];
            end
            c = c+length(yv);
        else
            temp_vect = f_HRTF_wrap(temp_vect);           
            m = 0;
           
            for new_col = 1:length(yv)
                m = m+1;
                Binaural_vectors{c_row,new_col}=[Binaural_vectors{c_row,new_col}
                    temp_vect(m,:)];
                Monaural_vectors{c_row,new_col}=[Monaural_vectors{c_row,new_col}
                    temp_mon(m,:)];
            end
        end % azn loop

    end
    %     =====================subtraction
    get_avg_binaural = [];
    for irow=1:size(Binaural_vectors,1)
        for icol = 1:size(Binaural_vectors,2)
            get_avg_binaural=[get_avg_binaural;Binaural_vectors{irow,icol}(end,:)];
        end
    end
    for irow=1:size(Binaural_vectors,1)
        for icol = 1:size(Binaural_vectors,2)
            Binaural_vectors{irow,icol}(end,:)=Binaural_vectors{irow,icol}(end,:)-mean(get_avg_binaural,1);
        end
    end
    if One1_Two2==2
        %     =====================subtraction
        get_avg_monaural = [];
        for irow=1:size(Monaural_vectors,1)
            for icol = 1:size(Monaural_vectors,2)
                get_avg_monaural=[get_avg_monaural;Monaural_vectors{irow,icol}(end,:)];
            end
        end
        for irow=1:size(Monaural_vectors,1)
            for icol = 1:size(Monaural_vectors,2)
                Monaural_vectors{irow,icol}(end,:)=Monaural_vectors{irow,icol}(end,:)-mean(get_avg_monaural,1);
            end
        end
    end
end % isub
% figure;semilogx(fsig(N1:N2),w(1:end-1),'.k-');
% xlim([f1 f2])
'enter'

if One1_Two2 == 1
    HRTF_cues{1} = Binaural_vectors;
else
    HRTF_cues{1} = Binaural_vectors;
    HRTF_cues{2} = Monaural_vectors;
end

out1 = f_LOOCV(One1_Two2,HRTF_cues,length(SN))
CM = out1.CM;
% for plotting:
figure
contourf(CM,'edgecolor','none')
m=colormap(gray);
colormap(flipud(m));

set(gca,'xtick',0:10:120,'ytick',0:10:120)
axis square
grid on
xlabel('Classified Speaker Index');
ylabel('True Speaker Index');
hold on;
% for i=10:10:120
%         plot([i i],[0 i],'k:');
%         plot([0 i],[i,i],'k:');    
% end


figure
for i=1:120
    for j=1:120
        if CM(i,j)>0
            s=plot(j,i,'s');hold on;
            set(s,'color','none','markersize',4,'markerfacecolor',(max(max(CM))+1-CM(i,j))./max(max(CM)).*[1 1 1])
        end
    end
end
title('Confusion Matrix', 'FontSize', 14)
set(gca,'xtick',0:10:120,'ytick',0:10:120)
xlabel('Classified Speaker Index');
ylabel('True Speaker Index');
axis square
grid on
exportgraphics(gcf, 'Figure_2_CM.png', 'Resolution', 300,'ContentType', 'vector')